 clear;
load('datastaticG')

phiI0=phiI;
phiX0=phiX;

for j=1:length(ss)-1
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);

zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;

ix=ixx(j);
x=xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);   
end
errorode(j+1)=0;


for j=flagstar-1:-1:1
%bs(j)=bs(j+1)-dbs(j+1)*ds;
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;

ix=ixx(j);
x=xx(j);

for kkj=1:600
db=(bs(j+1)-bs(j))/ds;

for kj=1:10
    
ixnew=(1-1/(bs(j)^(1-1/psi)*(A-ix-x)^(1/psi)/rho))/etaI;
ix=99/100*ix+ixnew/100;


end


errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);
while abs(errorode(j))>10^(-11)
derrorode=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi-1)*(1-1/psi)*(A-ix-x)/bs(j)^2*(-1))+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*((-1)/ds/bs(j)+(bs(j+1)-bs(j))/ds/bs(j)^2*(-1))+zetan(j)*(bsB(j)/bs(j))^(-gamma)*bsB(j)/bs(j)^2*(-1);
bs(j)=bs(j)-errorode(j)/derrorode;
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);   
end

%for kj=1:20
 xnew=(1-bs(j)/db*(1-etaI*ix))/etaX*s;
 x=9999/10000*x+xnew/10000;




end









ixx(j)=ix;
xx(j)=x;


end







for j=1:flagstar-1
s=ss(j);
betas=beta0+beta1*s;
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);
ix=ixx(j);
x=xx(j);    
    
db=(bs(j+1)-bs(j))/ds;
errori(j)=(1-1/(bs(j)^(1-1/psi)*(A-ixx(j)-xx(j))^(1/psi)/rho))/etaI-ixx(j);
%errorx(j)=(1-(bs(j)-s*db)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorx(j)=(1-bs(j)/db*(1-etaI*ixx(j)))/etaX*ss(j)-xx(j);
errorode(j)=rho/(1-1/psi)*(((A-ix-x)/bs(j))^(1-1/psi)-1)+(ix-etaI/2*ix^2-deltaI)-gamma*sigma^2/2+lambdas/(1-gamma)*(betas/(betas+1-gamma)-1)+((x/s-etaX/2*x^2/s^2-deltaX)-(ix-etaI/2*ix^2-deltaI))*s*(bs(j+1)-bs(j))/ds/bs(j)+zetan(j)/(1-gamma)*((bsB(j)/bs(j))^(1-gamma)-1);    
end
errori(j+1)=0;
errorx(j+1)=0;






continuedynamicG;
%generatemk;


% 
% 
% 
% subplot(3,2,1)
% plot(ss,bs,'LineWidth',1.5)
% hold on;
% 
% % plot(ss,bs0,'k:','LineWidth',1.5)
% % hold on;
% 
% plot(ss,bx,'r--','LineWidth',1.5)
% hold on;
% 
% 
% 
% 
% subplot(3,2,2)
% plot(ss,ixx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,ixx0,'k:','LineWidth',1.5)
% hold on;
% 
% subplot(3,2,3)
% plot(ss,xx,'LineWidth',1.5)
% hold on;
% 
% plot(ss,xx0,'k:','LineWidth',1.5)
% hold on;
% 
% 
% subplot(3,2,4)
% plot(ss(1:flagstar),errorode(1:flagstar),'LineWidth',1.5)
% hold on;
% plot(ss(1:flagstar),errori(1:flagstar),'r--','LineWidth',1.5)
% hold on;
% plot(ss(1:flagstar),errorx(1:flagstar),'k:','LineWidth',1.5)
% hold on;
% 
% subplot(3,2,5)
% plot(ss,ixx-ixx0,'LineWidth',1.5)
% hold on;
% 
% subplot(3,2,6)
% plot(ss,bs-bs0,'LineWidth',1.5)
% hold on;
% 
% % subplot(3,2,5)
% % plot(ss,dbs-dbsan,'LineWidth',1.5)
% % hold on;
